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ABSTRACT 


The  pseudospectral  method  is  especially  valuable  for  seismic  modeling  because  of  its  high  accuracy 
compared  to  other  numerical  techniques.  The  method  can  be  regarded  as  a  hmit  of  finite  difference 
of  increasing  orders,  and  a  process  of  trigonometric  interpolation,  thus  it  exhibits  high  accuracy. 
Stability  of  the  method  is  also  favorable.  Fourier  polynomials  are  especially  efficient  but  have  the 
disadvantage  of  forcing  periodicity,  and  Chebyshev  polynomials  are  somewhat  less  efficient  but 
are  more  flexible  in  application  of  boundary  conditions.  We  have  used  a  Fourier  pseudospectral 
method  in  the  horizontal  direction  and  Chebychev  polynomials  in  the  vertical  direction.  Curved 
grids  conforming  to  the  surface  topography  and  major  interfaces  are  made  possible  by  coordinate 
transformations.  A  full  viscoelastic  formulation  permits  convenient  implementation  of  attenuating 
layers  to  reduce  wrap-around  in  the  horizontal  direction.  The  result  is  an  efficient  method  for  2- 
and  3-D  linear  viscoelastic  wave  propagation. 
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1.  Objective 

1.1.  Introduction 

Anelasticity  and  anisotropy  of  the  earth  play  important  roles  in  wave  propagation,  especially  when  waves  travel  long 
distances.  Finite  difference  and  finite  element  methods  suffer  firom  inaccuracies  in  the  form  of  numerical  dispersion 
and  attenuation  that  make  it  difficult  to  simulate  wave  propagation  over  hundreds  to  thousands  of  wavelengths.  The 
pseudospectral  method  exhibits  very  striking  features  of  very  low  grid  density  and  high  efficiency,  compared  to  the 
standard  Cartesian  discrete  methods.  The  goal  here  is  to  simulate  wave  propagation  in  a  region  of  about  2,000  km 
in  azimuth  and  100-200  km  in  depth;  source  firequencies  are  0.5-10  Hz.  Since  wavelengths  range  from  120-4,000 
meters,  this  represents  propagation  of  about  500-20,000  wavelengths.  Most  FD  or  FE  schemes  either  exhibit  excessive 
numerical  artifacts  or  are  very  inefficient  for  such  problems. 

Due  to  space  limitations,  we  describe  only  the  isotropic  2-D  formulation.  However,  the  anisotropic  3-D  case  is  a 
straightforward  extension  and  a  detailed  analysis  is  available  firom  the  authors. 


1.2.  2-D  Viscoelastic  Wave  Equations 


When  we  talk  about  viscoelasticity,  there  are  mainly  two  different  ways  to  describe  the  strain-stress  relationship, 
i.e.  Voigt’s  model  or  Maxwell’s  model,  each  corresponding  to  connections  of  the  elastic  and  viscous  behavior  in 
parallel  and  series.  For  computational  consideration,  Voigt’s  model  is  highly  recommended  because  of  its  simplicity 
and  generality. 

For  an  isotropic-elastic  medium  the  stress-strain  relation  is 

an  =  XA+lJlcn,  i  =  1,2, 3,  (1) 

=  'pSij ,  J ,  J  =  1 , 2, 3,  (2) 

where  A  =  A  +  A'^  and  ]!=  fi  +  fi' -§^.11  is  beyond  the  scope  of  this  paper  to  include  a  derivation,  but  we  find  that  A' 
and  fi'  are  related  to  the  Q  for  P-  and  S-wave  propagation  by 


yJ{X  +  +  j(r2(A'  +  ln'fvj  +x+2^l 

2{X'  +  2fi')kvp 


(3) 


2fi'kv, 

We  see  that  in  this  model,  Qp  and  Q,  are  not  constant  but  are  functions  of  wavenumber.  The  viscoelastic  equations  of 
motion  are 

P^  =  ^[(A  +  ?r)^]  +  V-(7rVui)-H/i,  i  =  l,2,3,  (5) 

where  f  =  (/i ,  fz,  hV  is  the  applied  force. 

If  in  a  region  the  medium  is  partially  homogeneous  (/j  and  n'  are  constants),  then  in  this  region  we  have 


d'^e 


iX  +  2ji)Ae, 


(6) 


and 

where  ^  =  V  ■  u  and  w  =  V  x  u.  Thus  we  see  that  even  in  viscoelastic  media  the  dilation  wave  equation  is  similar  to  a 
“pure”  P-wave  equation  and  the  rotational  wave  equations  coincides  in  form  with  the  “pure”  S-wave  equation. 
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1.3.  2-D  Isotropic  Viscoelastic  Wave  Equations 


For  a  2-D  isotropic  viscoelastic  medium,  we  have 
d'^ui  d 


and 


d^U3 


dxi  [ 

d  ■ 
dx3 


,  _ , ,  <9ui  du3 


_  ,  /  dui  ,  du3 


_dui 


d  f_dui 


d  (_du3\  d  /_5u3^  r 


We  then  rewrite  these  equations  as  a  system  of  four  first  order  equations: 

dui 


dt 

dti3 

dt 


=  Vl, 


=  V3, 


(8) 

(9) 

(10) 

(11) 
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dvi  dv3 
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dxi  J  dx3  \  dx3  J 


(13) 


The  2-D  forward  modeling  problem  consists  of  solving  equations  (24-27)  in  rectangular  domains  with  a  free 
surface  on  X3  =  X5(xi).  The  equations  for  zero  normal  and  tangential  surface  stress  are  respectively: 


\5x3  dxi 

We  choose  radiation  boundary  conditions  at  the  bottom: 


(i+2^)gi+A|^  +  (A'*2p')^+A'gl-0, 


n  n  [«.(»■ .  *3)  3-  “3  a-cx. ,  + sin  aj(., ,  X,)  A  X  n  I 


Ui 

U3 


=  0. 


(14) 

(15) 

(16) 


where  |<?ij  |  <  7r/2,  j  =  1, 2,  •  ■  • ,  J.  The  radiation  boundary  conditions  (30)  are  imposed  such  that  the  wave  motions 
from  the  interior  of  the  domain  to  pass  through  the  boundary  X3  =  Xb{xi)  without  being  reflected.  In  general  it  is  not 
possible  to  find  practical  boundary  conditions  that  do  the  above  task  perfectly,  however,  the  artificial  reflections  can  be 
substantially  reduced.  Any  linear  combination  of  plane  waves  traveling  out  of  the  boundary  X3  =  Xb(xi)  at  angles  of 
incidence  ±1^1 ,  ±^2,  ■  •  • ,  with  speed  Va  satisfies  these  boundary  conditions  (Higdon,  1987),  where  the  absorption 
coefficients  (Meng  et  al.,  1989)  are 


^  iX+2py,  u’^y/pjX'  +  2p') 
2{X  +  2p)  2(A-f2^)3/2  ’ 


(17) 


_  p',  . 

‘  2p  2py^  ’ 


(18) 


No  explicit  boundary  conditions  are  use  on  the  sides.  The  natural  (periodic)  conditions  are  employed  for  Fourier 
polynomials.  Absorption  is  accomplished  by  using  hin  absorbing  layers  with  p'  and  X'  gradually  increasing  towards 
the  boundaries. 
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To  obtain  a  rectangular  geometry  for  a  problem  with  an  irregular  top  and  bottom  (due  to  surface  topography  and 
curvature  of  the  earth)  we  employ  a  transformadon  from  (ii ,  X3)  to  (xi,  77),  where  77  =  77(xi ,  X3)  (Fig.  1).  The  equations 
of  motion  then  become 


+  Vx,- 


y  +  u' 


dvi  dvi  dv-i 

dxi  dr)  dr] 


(  d  5  '\  r  ,  /  (9u3  5t;3  \  1  f  d\\,f  dv^\ 

The  free  surface  boundary  conditions  at  the  top  become: 

o  \{  A dux\  ,  f  dv3\  .,fdvi  t 

(A  +  27i)(^77.3-^j+A(^— +  77..— j+(A  +2/z)(^77.,— j+A 

/  5ui  du3  du3\  ,{  dv\  dv3  dv3\  ^ 

The  radiation  boundary  conditions  at  the  bottom  take  the  form 


(22) 

(23) 

(24) 


/ dui  dui  \  .  (  du\  \  vi  ^ 


(25) 


Q:(Xi,X3)U3  +  COS0(Xi,X3)(  — +  77,,—  1  +  Sin  (^(Xi,  X3)  I  T?,, — 


+  T^  =  0. 
i;(xi,X3) 


These  equations  are  to  be  approximated  through  collocation  and  quadrature  (we  use  fourth  order  Runge-Kutta). 
The  solutions  are  represented  as  2-D  discrete  Fourier-Chebyshev  transforms  in  (xi,  77).  Differentiations  with  respect 
to  xi  and  77  are  carried  out  in  the  usual  way.  The  77  derivatives  of  m,  U3,  v\  and  vi  result  in  four  unknovras  per 
value:  the  highest  Chebyshev  modes  are  lost  in  the  differentiation  process.  In  principal,  these  are  to  be  determined 

by  the  47V3:  top  and  bottom  boundary  conditions.  Due  to  the  higher  derivatives  in  the  equations  of  motion,  there  are 
actually  a  total  of  lANx  unknown  highest  Chebyshev  modes.  The  highest  modes  in  the  equations  are  irrelevant  and  are 
ignored.  Instead,  the  highest  modes  in  the  AN^  imknowns  are  updated  according  to  the  ANx  top  and  bottom  boundary 
conditions. 

Stability  of  Runge-Kutta  for  this  problem  has  been  studied.  Details  are  available  from  the  authors. 
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1.4.  About  the  Grids 

The  strategy  adopted  for  choosing  the  mapping  function  77  is  to  flatten  the  surface  and  major  contours.  “Major  contours” 
are  those  interfaces,  typically  the  basement  and  moho,  that  are  continuous  across  the  section.  Fig.  I  shows  one  of  a 
set  of  “deep  sedimentary  basin”  models  that  are  being  tested.  Fig.  2  shows  a  model,  derived  from  seismic  data,  for  a 
trans-Urals  DSS  path. 

Once  the  “major  contours”  are  identified,  the  average  travel  time  between  the  contours  is  computed  (p-waves 
have  been  used  so  far),  and  constant  values  of  rj  are  assigned  to  the  major  contours  so  that  the  differences  in  p  are 
proportional  to  the  travel  time.  In  this  way,  increments  of  constant  At)  correspond  roughly  to  equal  increments  of  travel 
time.  This  is  done  so  as  to  improve  the  stability  properties  of  the  time  integration. 


1.5.  Preliminary  Results 

Preliminary  calculations  indicate  that  it  is  practical  to  simulate  wave  propagation  for  2-D  models  of  about  2,000  by 
200  km.  Such  simulations  require  several  days  of  compuation  on  a  SGI  Power  Challenge  (300  MFlops). 

We  are  in  the  process  of  quantifying  the  effects  of  particular  types  of  crustal  features  on  regional  waveforms. 
Features  under  study  include  deep  sedimentary  basins  (see  Fig.  1),  moho  uplift  and  lateral  variation  in  surface  features. 
Also,  some  realistic  models  (such  as  the  Soviet  DSS  model  in  Fig.  2)  are  currently  being  tested. 


2.  Recommeadations  and  Future  Plans 

Other  authors  have  performed  simulations  for  particular  models  (e.g.,  Jih  (1 994)  uses  a  finite  difference  method);  there 
is  not  yet  completed  a  comprehensive  and  quantitative  study  of  the  effects  of  crustal  features.  Our  plan  is  to  use  a 
series  of  models,  for  example  the  deep  sedimentary  basin  models,  varying  the  size  and  elastic  parameters  of  the  basin, 
to  determine  how  the  effects  of  such  basins  on  crustal  wave  propagation,  such  as  attenuation  of  L^,  depend  on  the 
parameters  of  the  basin.  Additionally,  we  will  study  the  effects  of  moho  uplift,  combined  basin  and  uplift,  and  surface 
topographical  roughness  and  other  surface  features. 

One  weakness  of  this  method  is  that,  because  of  its  2-D  nature,  amplitudes  and  decay  times  are  not  correct,  and 
relative  amplitudes  of  differing  phases  are  not  readily  obtainable.  Although  the  method  described  in  this  paper  can 
readily  be  applied  to  the  3-D  problem,  it  is  apparent  that  elastic  wavefield  computations  in  a  region  of  2,000  by  2,000 
by  200  km  are  not  practical  with  any  computer  available  to  the  authors.  However,  two  acceleration  techniques  may  be 
employed  when  solving  the  problem  of  wave  propagation  between  two  points  (source  and  receiver): 

1.  Model  only  a  narrow  channel  between  the  source  and  receiver  (see  Fig.  3).  Use  absorbing  boundary  conditions  on 
the  lateral  boundaries  to  avoid  reflections  or  wrap-around. 

2.  Use  smaller  tracking  grids  that  follows  the  major  wave  fronts.  A  full  simulation  would  then  require  several  runs, 
each  to  capture  the  modes  with  essentially  different  paths,  but  each  run  would  be  much  quicker  than  a  global 
computation. 
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Fig.  1.  Deep  sedimentary  basin  model,  showing  eta  grid  lines 
conforming  to  the  major  contours,  and  eta  spacing  proportional 
to  travel  time. 
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Fig.  2.  Trans-Urals  Soviet  DSS  model. 
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